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We explore the newly discovered "hangup-kick" effect, which greatly amplifies the recoil for con- 
figurations with partial spin-/ orbital- angular momentum alignment, by studying a set of 48 new 
simulations of equal-mass, spinning black-hole binaries. We propose a phenomenological model for 
the recoil that takes this new effect into account and then use this model, in conjunction with statis- 
tical distributions for the spin magnitude and orientations, based on accretion simulations, to find 
the probabilities for observing recoils of several thousand km s~ . In addition, we provide initial 
parameters, eccentricities, radiated linear and angular momentum, precession rates and remnant 
mass, spin, and recoils for all 48 configurations. Our results indicate that surveys exploring peculiar 
(redshifted or blueshifted) differential line-of-sight velocities should observe at least one case above 
2000 km s _1 out of four thousand merged galaxies. On the other hand, the probability that a 
remnant BH recoils in any direction at a velocity exceeding the ~ 2000 km s _1 escape velocity of 
large elliptical galaxies is 0.03%. Probabilities of recoils exceeding the escape velocity quickly rise 
to 5% for galaxies with escape velocities of 1000 km s _1 and nearly 20% for galaxies with escape 
velocities of 500 km s _1 . In addition the direction of these large recoils is strongly peaked toward 
the angular momentum axis, with very low probabilities of recoils exceeding 350 km s _1 for angles 
larger than 45° with respect to the orbital angular momentum axis. 

PACS numbers: 04.25.dg, 04.30.Db, 04.25.Nx, 04.70.Bw 



I. INTRODUCTION 

Speculations about the relevance of gravitational re- 
coils in astrophysical black-hole binary (BHB) mergers 
can be traced back at least thirty years [TJ |2] . The crucial 
scale of the problem is when those recoils reach veloci- 
ties comparable to the escape velocities of the relevant 
structures, i.e. globular clusters, which have escape ve- 
locities of 10s of km s _1 , and dwarf, spiral, and giant 
elliptical galaxies, which have escape velocities from 100s 
to ~ 1000 km s _1 for normal galaxies. For large galaxies 
undergoing major mergers the effective escape velocity 
can be up to a factor of a few higher at the time of coa- 
lescence, as the central potential well deepens rapidly at 
that time. Once the merger is complete and the stellar 
systems begins to relax, the potential becomes shallower 

Early attempts to compute recoil velocities from BHB 
mergers used perturbative [H [5] and post-Newtonian ap- 
proximations (see [6] for a review up to 2005) and found 
recoils up to a few hundred km s _1 , but uncertainties in 
the computations were of the same order of magnitude as 
those velocities (see [7] for a more current review). The 
first computation that used full numerical simulations 
within the Lazarus approach produced similar results [5] . 

The accurate computation of recoil velocities had to 
wait for the 2005 breakthroughs [5HTT] in Numerical rel- 
ativity (NR), since it proved to be a genuinely strong- 
field, highly nonlinear General Relativistic phenomenon. 
The first systematic study of recoil velocities considered 



unequal mass, nonspinning BHBs 12J. That study found 
that the maximum recoil velocity for non-spinning BHBs 
is 175 km s , which occurs for a mass ratio near 1:3. 

Unexpectedly, spinning BHBs, with individual spins 
anti-aligned with each other and both parallel to the an- 
gular momentum direction were found to produce recoils 
of over a factor two larger than the unequal-mass maxi- 
mum [T31 [Uj , and a revolution occurred when it was dis- 
covered [15] that a configuration of the spins lying in the 
orbital plane led to recoils of almost [THIE] 4000 km s _1 . 
This last figure caught the attention of observational as- 
tronomers who began to look for these highly-recoiling 
BHs by searching the spectral data of galaxies for dif- 
ferential redshifts of several thousand km s _1 . The idea 
there was that gas close to the BH would remain bound to 
it, while gas further out in the accretion disk would be left 
behind. The two gas components would then have differ- 
ent relative redshifts. Initial searches produced the first 
supermassive recoiling BH candidates 18 -21] and now 
more thorough surveys have increased the numbers of 
potential candidates to several dozen [SU [53] . These ob- 
servations may provide the first confirmation of a general 
relativistic strong field, highly-dynamical, full-numerical 
prediction. 

A recent study |17j pointed out that configurations 
with partially aligned spins, which we call the "hangup- 
kick" configuration, can lead to even larger recoil veloc- 
ities, of nearly 5000 km s _1 . More importantly, these 
configurations are favored with respect to the "spin in 
the orbital plane" configuration by the effects of accre- 
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tion on the BHs during very early orbital (Newtonian) 
stages [24j [25] . We address this question in more detail 
in this paper. 

This paper is organized as follows. In Sec. [IT] we review 
the numerical methodology to perform the simulations. 
In Sec. |III[ we describe the initial configurations of a fam- 
ily of BHBs chosen to model "hangup-kicks". In Sec. |IV| 
we provide the main results of these evolutions in ta- 
bles of radiated energies, angular and linear momenta, 
as well as final remnant mass and spin. We then model 
recoils using empirical fitting formulas. In Sec. [V] we de- 
scribe smoothed particle hydrodynamics (SPH) simula- 
tions that model accretion onto BHBs to obtain the spin 
magnitude and direction distribution of the individual 
spins in merging BHBs for our full numerical simulations. 
Using these spin direction and magnitude distributions, 
we give predictions for the recoil distribution and the 
probabilities of observing large recoils. We discuss the 
consequences and future extensions of these techniques 
in Sec. [VHP 

II. NUMERICAL RELATIVITY TECHNIQUES 

We use the TwoPunctures thorn [55] to generate initial 
puncture data [37] for the black-hole binary (BHB) simu- 
lations described below. These data are characterized by 
mass parameters m p , which are not the horizon masses, 
of each BH, as well as the momentum and spin of each 
BH. We evolve these BHB data-sets using the LazEv [25] 
implementation of the moving puncture approach |10[lllj 
with the conformal function W = ^/x = cxp(— 2</>) sug- 
gested by Ref. For the runs presented here, we use 
centered, eighth-order finite differencing in space |30j and 
a fourth-order Runge Kutta time integrator. (Note that 
we do not upwind the advection terms.) 

Our code uses the Cactus/EinsteinToolkit [3T1I32"] 
infrastructure. We use the Carpet [33] mesh refinement 
driver to provide a "moving boxes" style of mesh refine- 
ment. In this approach refined grids of fixed size are ar- 
ranged about the coordinate centers of both holes. The 
Carpet code then moves these fine grids about the com- 
putational domain by following the trajectories of the two 
BHs. 

We obtain accurate, convergent waveforms and horizon 
parameters by evolving this system in conjunction with a 
modified 1+log lapse and a modified Gamma-driver shift 
condition [101 (Ml I3S], and an initial lapse a(t = 0) = 
2/(l + tpj BL ), where ipBL is the Brill-Lindquist conformal 
factor and is given by 

71 

V> B L = i + X>i/( 2 lr-n|), 

where r, is the coordinate location of puncture i. The 
lapse and shift are evolved with 

(d t - li l di)a = -2aK, (la) 
d t (3 a = (3/4)f a - V (3 a , (lb) 



where we use r\ = 2 for all simulations presented below. 

We use AHFinderDirect [36] to locate apparent 
horizons. We measure the magnitude of the horizon spin 
using the Isolated Horizon algorithm detailed in Ref. [37] . 
Note that once we have the horizon spin, we can calculate 
the horizon mass via the Christodoulou formula 

m H = y< + Sy(K r ), (2) 

where m irr = \J Aj (167r) and A is the surface area of 
the horizon, and Sh is the spin angular momentum of 
the BH (in units of M 2 ). In the tables below, we use 
the variation in the measured horizon irreducible mass 
and spin during the simulation as a measure of the error 
in these quantities. We measure radiated energy, linear 
momentum, and angular momentum, in terms of the ra- 
diative Weyl Scalar ^4, using the formulas provided in 
Refs. [35] I2S]- However, rather than using the full ^4, 
we decompose it into i and m modes and solve for the 
radiated linear momentum, dropping terms with I > 5. 
The formulas in Refs. [351 SS] are valid at r = 00. We ex- 
tract the radiated energy-momentum at finite radius and 
extrapolate to r = 00 using both linear and quadratic 
extrapolations. We use the difference of these two ex- 
trapolations as a measure of the error. 

III. SIMULATIONS 

We evolved a set of 48 equal-mass, spinning, qua- 
sicircular, "hangup-kick" configurations, with 30 simu- 
lations having individual BH spins of magnitude a = 
1/V2 and 18 simulations having BH spin magnitudes of 
a = 0.9, where a is the dimensionless spin of the BH 
(a = Sh/Mjj, where Sh is the spin angular momentum 
and Mh is the mass of the BH). In the "hangup- kick" 
configuration, the z components of the individual spins 
are equal, while the projections of the individual spins 
onto the orbital plane are equal in magnitude but oppo- 
site in direction. The a = 1/ V2 configurations were split 
into five sets of 6, where the runs in each individual set 
had the same initial angle 9 between the spin direction 
and orbital angular momentum direction (here we chose 
9 = 22.5°, 45°, 60°, 120°, 135°). In each set with as 
given 9, we chose the initial orientation between the 
in-plane spin and linear momentum to be 0°, 30°, 90°, 
130°, 210°, and 315°. For the a = 0.9 runs, we used 
the same initial 6 cf>i configurations for 9 — 60°, 9 = 30°, 
and 9 = 15°. We combine these results with the simu- 
lations of [50J (which have 9 = 90°) in order to perform 
our analysis below. 

Initial data parameters for the 48 simulations are 
given in Table [I] We denote these configurations by 
AsTHxxxPHyyy, where s indicates the approximate in- 
dividual spin magnitude (7 for en = l/-\/2 and 9 for 
oti = 0.9), xxx indicates the angle the spin makes with 
respect to the z axis, and yyy indicates the angle the 
spin makes with respect to the y axis. Here xxx and yyy 
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FIG. 1: xy plane projections of the trajectories for various 
"hangup-kick" configurations. (Top Left) Trajectory for the 
A7TH22.5PH0 configuration, (Bottom Left) trajectory for the 
A7TH135PH0 configuration, (Top Right) trajectory for the 
A9TH15PH0 configuration, (Bottom Right) trajectory for the 
A9TH60PH0 configuration. The plot shows the trajectories 
for configurations with the largest and smallest inclination 
angle for the a = 1/V2 and a = 0.9 configurations. Note that 
the eccentricity is larger for large 9 and that the eccentricity 
decreases more slowly. 



are in degrees. An eccentricity reduction procedure like 
those given in J4TJ 02] could be used to generate config- 
urations with very low eccentricity, but the amount of 
time required to reduce the eccentricity for 48 configu- 
rations would have been too long. We therefore chose 
to start from 3.5 Post-Newtonian (PN) quasicircular or- 
bital parameters from further separations, such that each 
binary completed 5-6 orbits prior to merger, and then re- 
lied on the radiation of angular momentum during this 
5-6 orbit inspiral to reduce the eccentricity. The initial 
separations varied between 10.16M and 8.2M, depend- 
ing on the magnitude of the hangup effect, with smaller 
initial separations for configurations that exhibit larger 
hangups. Example trajectories for several of these con- 
figurations are given in Figs. [T] and [2j 



The orbital motion of these BHBs has an interesting 
property, the spin precession frequency and the orbital 
frequency agree right near merger (in these coordinates). 
Initially the spin precession frequencies are much lower 
than the orbital frequency, but ramp up dramatically 
near merger. Figure [3] shows this behavior for the dif- 
ferent 9 and a configurations (for clarity in the plot, we 
only show the <pi — configurations). Interestingly, the 
strongest effect seems to be due to the inclination angle 
9 (which is measured with respect to the orbital angu- 
lar momentum axis, i.e. the z axis) rather than on the 
projected z spin or total spin. For a given a and 9 there 
are variations with (/>,;, but these are smaller than the 
variations with 6. This rapid increase of the precession 
frequency near merger (as measured with the techniques 
of |43] ) . increasing all the way up to the orbital frequency, 
is in contrast with the much milder increase in the spin 
magnitude due to weak tidal effects [33]. It also lends sup- 
port to modeling of black hole merger assuming geodetic 
precession [43] , 

In Table [Tl] we compare the radiated mass and angular 
momentum as calculated directly from ip4 to the corre- 
sponding quantities derived from the remnant mass and 
spin. The difference between these quantities is a bet- 
ter measurement of the error than those derived from 
variations in the final horizon mass and spin or in ex- 



trapolation of the waveform to infinity. That is, there 
are systematic errors due to truncation error and finite 
extraction radius, and the difference between these two 
measurements gives a lower bound to the error. 

The dimensionless spin a mcrge r, the orientation tp of 
the spin during the final orbit and plunge, as well as 
the remnant BH properties, including recoil velocity, are 



given in Table III Note here that the orientation tp is 
the angle that the spin of BH1 (the BH originally lo- 
cated on the positive x axis) makes with the spin of 
BH1 in the corresponding AsTHxxxPHO configuration in 
a rotated frame where infall directions all coincide (see 
[46] ) . Also note that the largest measured recoil for these 
runs is (4079.5 ± 10.1) km s" 1 for the A9TH60PH30 
and A9TH60PH210 configurations, which exceeds both 
the largest measured quasicircular "superkick recoil" of 
3300 km s _1 [47] and even the theoretical "superkick 
maximum" recoil of 3680 ± 130 km s _1 [30]. The values 
given for the spins near merger should only be taken as 
an approximation. The spin-up apparent in the A9TH60 
simulations (i.e. the difference between a me rger and the 
initial spin of a — 0.9) was due to the lower resolution 
used in these simulations (simulations of A9TH45 with 
the same resolution as A9TH60 showed an even stronger 
spin-up at late times that converged away with higher 
resolution). Interestingly, these highly-spinning configu- 
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TABLE I: Initial data parameters for the 48 "hangup-kick" configurations. In all cases the puncture masses were chosen such 
that the total ADM mass of the binary was 1.0 ± 1CP 6 M. Here the punctures are located at ±{x, 0, 0) with momenta ±(0,p, 0) 
and spins S = (±S x ,±S y , S z ). The approximate initial eccentricities, eccentricities measured over the last orbit, and the 
number of orbits, are also given. 
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rations can exhibit both spin-up and spin-down, depend- 
ing on the location of the refinement boundaries, when 
not fully resolved. The other A9 runs used higher reso- 
lution and show much better spin conservation. 

We evolve the A9TH15 and A9TH30 configurations 
with 10 levels of refinement and maximum resolution of 



h = M/153.6. The width of this level was 2 x 0.35M, 
while the radius of the horizons grew to 0.24M. Our 
initial explorations used grids that were smaller in ra- 
dius, but we found that using larger grids improved the 
spin conservation considerably. The A9TH60 configura- 
tion were evolved with grids a factor of 1.2 coarser, and 
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consequently the spins near merger are not as accurate. 
Note that these A9TH60 runs were performed first, and 
based on the errors in these simulations, we refined the 
grid for the other A9 runs. 



FIG. 2: The elevation of the trajectory as a function of time 
for the several of the "hangup-kick" configurations. Note that 
the "bobbing" amplitude does not necessarily correspond to 
a large recoil. The A7TH135PH0 configuration has a factor 
of 2 smaller recoil than he A7TH60PH0 configuration, but a 
slightly larger bobbing amplitude. 



IV. RESULTS AND MODELING OF RECOIL 
VELOCITIES 



With the discovery of very large recoil [15j velocities 
for certain configurations of merging spinning BHBs, the 
need for an empirical model for the recoil velocity as 
a function of the progenitor's parameters was apparent. 
Our approach to provide that phenomenological formula 
was based on the observation that the recoil of spinning 
BHs is largely generated around the time of merger of the 
two holes [H]; and that this nearly instantaneous burst 
of radiation of linear momentum can be modeled by a 
parametrized dependence of the leading (on spins and 
mass ratio) post-Newtonian (PN) expressions for the lin- 
ear momentum radiated 149 . 



In Ref. |50j we extended our original empirical formula 
for the recoil velocity imparted to the remnant of a BHB 
merger [TS1 [TB] to include next-to-leading-order correc- 
tions (based on the PN work of [H]), still linear in the 
spins 



where 
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and 77 = <?/(l + q) 2 , with q = mi/7712 the mass ra- 
tio of the smaller to larger mass hole, cJj = Si/m 2 , 
rrii is shorthand for mm the mass of BH i, the index 
_L and || refer to perpendicular and parallel to the or- 
bital angular momentum respectively, ei,e2 are orthog- 
onal unit vectors in the orbital plane, and £ measures 
the angle between the unequal mass and spin contribu- 
tion to the recoil velocity in the orbital plane. The an- 
gles 0a and 4>s are defined as the angle between the 
in-plane component A 1 - = M(S% jm<i — Si/m\) and 
S 1 - = Sf- + S% respectively and a fiducial direction at 
merger (see Ref. [35] for a description of the technique). 
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TABLE II: A comparison of the radiated mass and angular momentum with the predictions based on the final remnant mass 
and spin and the initial ADM mass and angular momentum for the 48 "hangup-kick" configurations. 



CONF 


(Madm - 


- M H )/M 


8E TaA /M 






(Vadm 


S 


h)/M 2 


8JLJM 


2 




A7TH22.5PH0 


0.065949 


-1- 




0.063479 


± 


0.000169 


0.452347 


± 


0.001590 


0.440182 


± 


0.006183 


A7TH22.5PH30 


0.065440 


4- 

ZL 


n nnm 9^ 


0.063034 


zb 


0.000164 


0.450602 


zb 


0.001706 


0.438756 


zb 


0.006006 


A7TH22.5PH90 


0.065188 


± 


0.000121 


0.062822 


zb 


0.000171 


0.450332 


zb 


0.001675 


0.439398 


zb 


0.005161 


A7TH22.5PHI30 


0.065842 


-|- 


o nnm n5 

u . uuu J. uu 


0.063363 


± 


0.000190 


0.452581 


zb 


0.001493 


0.440443 


zb 


0.005901 


A7TH22.5PH210 


0.065440 


-|- 


n nnm 25 

u • uuu lzjO 


0.063035 


± 


0.000165 


0.450601 


zb 


0.001706 


0.438766 


± 


0.005996 


A7TH22.5PH3I5 


0.065910 


-|- 


n nnm n5 

U . UUU 1UU 


0.063423 


± 


0.000188 


0.452761 


± 


0.001487 


0.440460 


± 


0.006077 


A7TH45PH0 


0.058764 


4- 

ZL 


n ooooofi 


0.056538 


± 


0.000181 


0.417697 


± 


0.000084 


0.403669 


± 


0.007395 


A7TH45PH30 


0.058594 


-|- 


n ooooos 

U.UUUUUO 


0.056460 


zb 


0.000153 


0.416460 


zb 


0.000088 


0.403290 


zb 


0.006741 


A7TH45PH90 


0.056266 


± 


0.000010 


0.054394 


± 


0.000120 


0.408936 


zb 


0.000072 


0.397875 


zb 


0.004794 


A7TH45PH130 


0.056959 




o nnnm n 

U • UUUU -L U 


0.054940 


zb 


0.000155 


0.412016 


zb 


0.000078 


0.399658 


zb 


0.005982 


A7TH45PH210 


0.058594 




n nnnnns 

U . UUUUUO 


0.056460 


zb 


0.000153 


0.416459 


zb 


0.000088 


0.403290 


zb 


0.006739 


A7TH45PH3I5 


0.057167 


± 


n nnnoos 


0.055120 


zb 


0.000159 


0.412756 


zb 


0.000078 


0.400125 


± 


0.006224 


A7TH60PH0 


0.050887 


4- 


u.uuuuuu 


0.049046 


zb 


0.000145 


0.379033 


± 


0.000032 


0.365962 


± 


0.007685 


A7TH60PH30 


0.052125 


-|- 


00000^ 

U.UUUUUO 


0.050267 


zb 


0.000146 


0.383455 


zb 


0.000029 


0.371173 


zb 


0.006948 


A7TH60PH90 


0.051234 


± 


0.000005 


0.049628 


zb 


0.000094 


0.378829 


zb 


0.000027 


0.368250 


zb 


0.005342 


A7TH60PH130 


0.049522 


± 


000004 


0.047918 


zb 


0.000103 


0.372896 


zb 


0.000025 


0.361386 


± 


0.006580 


A7TH60PH210 


0.052124 


-|- 


o oonnns 

u . uuuuuo 


0.050266 


zb 


0.000146 


0.383456 


zb 


0.000029 


0.371168 


± 


0.006945 


A7TH60PH315 


0.049445 


-|- 


o oonnn4 

U . UUUUU4 


0.047811 


zb 


0.000112 


0.372866 


± 


0.000030 


0.360951 


± 


0.006697 


A7THI20PH0 


0.033193 


4- 

ZL 


000009 

U.UUUUUZ 


0.032297 


zb 


0.000048 


0.303865 


± 


0.000005 


0.298884 


zb 


0.006258 


A7THI20PH30 


0.033107 


-|- 


000009 
u.uuuuuu 


0.032297 


zb 


0.000010 


0.303095 


zb 


0.000005 


0.295666 


± 


0.008881 


A7TH120PH90 


0.031686 


± 


0.000002 


0.031037 


zb 


0.000015 


0.297018 


zb 


0.000005 


0.292110 


± 


0.006177 


A7THI20PH130 


0.031974 


-|- 


o nnonn? 

U . UUUUUO 


0.031188 


zb 


0.000044 


0.298771 


zb 


0.000005 


0.297030 


± 


0.002584 


A7THI20PH2I0 


0.033107 


-|- 


o oonnn? 

U . UUUUUO 


0.032298 


± 


0.000010 


0.303097 


zb 


0.000005 


0.295663 


± 


0.008889 


A7THI20PH3I5 


0.032117 


-|- 


o ooonn? 

U . UUUUUZi 


0.031308 


± 


0.000051 


0.299419 


± 


0.000005 


0.297731 


± 


0.002485 


A7THI35PH0 


0.029675 


-|- 


000004 
u.uuuuuu 


0.029004 


zb 


0.000024 


0.287481 


zb 


0.000006 


0.290758 


zb 


0.000188 


A7THI35PH30 


0.030091 


4- 

ZIZ 


000004 
u.uuuuuu 


0.029449 


± 


0.000043 


0.289543 


zb 


0.000005 


0.293205 


± 


0.000038 


A7THI35PH90 


0.030000 


± 


0.000004 


0.029437 


± 


0.000063 


0.288696 


zb 


0.000006 


0.294013 


± 


0.001881 


A7THI35PH130 


0.029415 




o oonnns 

u . uuuuuo 


0.028803 


zb 


0.000032 


0.285893 


zb 


0.000006 


0.290818 


zb 


0.001668 


A7THI35PH2I0 


0.030085 


± 


000004 


0.029448 


zb 


0.000048 


0.289513 


zb 


0.000006 


0.293050 


zb 


0.000149 


A7THI35PH3I5 


0.029385 


-1- 


o oonnn4 

U . UUUUU4 


0.028766 


± 


0.000032 


0.285792 


± 


0.000006 


0.290244 


± 


0.001187 


A9THI5PH0 


0.086926 


-1- 


000499 

U . UUU4iZ 


0.082160 


± 


0.000612 


0.540966 


± 


0.004231 


0.523516 


± 


0.003980 


A9THI5PH30 


0.087312 


4- 
zn 


OOO^Qfi 


0.082564 


zb 


0.000612 


0.541701 


zb 


0.003972 


0.524385 


zb 


0.004228 


A9THI5PH90 


0.086773 


zb 


0.000341 


0.082188 


zb 


0.000602 


0.539095 


zb 


0.003287 


0.523929 


zb 


0.003074 


A9THI5PH130 


0.086317 


-|- 


000384 


0.081695 


zb 


0.000599 


0.538447 


zb 


0.003722 


0.522743 


zb 


0.002983 


A9THI5PH210 


0.087315 


-|- 


000394 


0.082564 


zb 


0.000612 


0.541724 


zb 


0.003950 


0.524384 


± 


0.004228 


A9THI5PH3I5 


0.086325 


-|- 


000394 

U . UUUOC/4: 


0.081686 


± 


0.000600 


0.538588 


± 


0.003835 


0.522588 


± 


0.003187 


A9TH30PH0 


0.077615 


-|- 


OOO^Q 

U.UUUOOc/ 


0.073045 


± 


0.000461 


0.509243 


± 


0.005934 


0.484306 


± 


0.004248 


A9TH30PH30 


0.076957 


-|- 


000^1 fi 

U.UUUOIU 


0.072527 


zb 


0.000493 


0.507180 


zb 


0.005340 


0.483676 


zb 


0.004906 


A9TH30PH90 


0.078900 


zb 


0.000166 


0.074900 


zb 


0.000485 


0.510059 


zb 


0.002727 


0.496097 


zb 


0.000678 


A9TH30PH130 


0.079891 


zb 


0.000209 


0.075617 


zb 


0.000486 


0.513344 


zb 


0.003666 


0.496170 


zb 


0.000939 


A9TH30PH2I0 


0.076956 


zb 


0.000316 


rx rx f r~ 

0.072556 


zb 


0.000469 


0.507180 


zb 


0.005341 


0.483678 


zb 


0.004902 


A9TH30PH3I5 


0.079810 


zb 


0.000224 


0.075477 


zb 


0.000486 


0.513380 


± 


0.003973 


0.495112 


zb 


0.001297 


A9TH60PH0 


0.056289 


± 


0.000311 


0.055634 


zb 


0.000322 


0.387215 


± 


0.007605 


0.411510 


± 


0.009728 


A9TH60PH30 


0.059774 


± 


0.000325 


0.058422 


zb 


0.000326 


0.402796 


zb 


0.005399 


0.422236 


zb 


0.007268 


A9TH60PH90 


0.054809 


zb 


0.000087 


0.052769 


zb 


0.000184 


0.394238 


zb 


0.001762 


0.396247 


zb 


0.003734 


A9TH60PH130 


0.055249 


± 


0.000168 


0.053567 


zb 


0.000241 


0.392036 


zb 


0.003341 


0.399533 


± 


0.008047 


A9TH60PH210 


0.059774 


± 


0.000325 


0.058422 


zb 


0.000326 


0.402796 


zb 


0.005399 


0.422236 


± 


0.007268 


A9TH60PH315 


0.054902 


± 


0.000201 


0.053558 


zb 


0.000260 


0.388072 


zb 


0.004485 


0.400454 


zb 


0.008515 



Note that A = M(5 l 2/m2 — Si /mi) can be expressed as 
A = M 2 (d?2 — qdi) Phases </>i and 4>2 depend on 

the initial separation of the holes for quasicircular orbits 
(astrophysically realistic evolutions of comparable masses 
BHs lead to nearly zero eccentricity mergers). 

The most recent published estimates for the above pa- 
rameters can be found in \46\ 52 and references therein. 



The current best estimates are: A m = 1.2 x 10 4 km s 1 , 
B m = -0.93, H = (6.9 ± 0.5) x 10 3 km s" 1 , K = 
(5.9±0.1)xl0 4 km s"\ and £ - 145°, and K s = -4.254. 
Here we set Bjj and Bk to zero, which is consistent with 
the findings in [SD], where it was found that the uncer- 
tainties in the coefficients are of the same magnitude as 
the coefficients themselves. 

Although the post-Newtonian approximation fails to 
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TABLE III: Merger and remnant BH properties of the 48 configurations. Sh is the spin angular momentum of the remnant, 
Mh is the Christodoulou mass, V r e CO ii is the recoil velocity, a m erger is an approximate value of the dimensionless spin during 
the merger phase, and is the angle between the direction of the spin of BH1 (in the rotated frame) and the spin of BH1 in 
the corresponding PHO configuration (see Section IV). 



CONF 


M H /M 


Sh/M 2 


K-o r oii(km s 1 ) 

1 (. (. Oil \ / 


^■merger 


ip 


A7TH22.5PH0 


0.934051 ± 0.000113 


0.760575 ± 0.001590 


-925.3 ± 1.0 


0.71 





A7TH22.5PH30 


0.934560 ± 0.000125 


0.762320 ± 0.001706 


-7.9 ±2.3 


0.71 


31.53 


A7TH22.5PH90 


0.934812 ± 0.000121 


0.762590 ± 0.001675 


1531.5 ±2.6 


0.71 


91.68 


A7TH22.5PH130 


0.934158 ± 0.000105 


0.760341 ± 0.001493 


1735.2 ± 1.2 


0.71 


131.14 


A7TH22.5PH210 


0.934560 ± 0.000125 


0.762321 ± 0.001706 


8.0 ± 2.3 


0.71 


211.52 


A7TH22.5PH315 


0.934089 ± 0.000105 


0.760161 ±0.001487 


-1703.7 ± 1.1 


0.71 


315.97 


A7TH45PH0 


0.941235 ± 0.000006 


0.730165 ± 0.000084 


2527.7 ± 5.4 


0.71 





A7TH45PH30 


0.941406 ± 0.000008 


0.731403 ± 0.000088 


1708.6 ± 1.2 


0.71 


29.33 


A7TH45PH90 


0.943734 ±0.000010 


0.738926 ± 0.000072 


-1199.6 ±4.8 


0.71 


91.22 


A7TH45PH130 


0.943041 ± 0.000010 


0.735846 ± 0.000078 


-2486.4 ±0.2 


0.71 


131.57 


A7TH45PH210 


0.941406 ± 0.000008 


0.731403 ± 0.000088 


-1708.2 ± 1.3 


0.71 


209.35 


A7TH45PH315 


0.942833 ± 0.000008 


0.735106 ± 0.000078 


2569.6 ±0.7 


0.71 


316.39 


A7TH60PH0 


0.949113 ± 0.000004 


0.706585 ± 0.000032 


-2786.0 ± 4.6 


0.71 





A7TH60PH30 


0.947875 ± 0.000003 


0.702163 ± 0.000029 


-2886.8 ±2.9 


0.71 


27.78 


A7TH60PH90 


0.948766 ± 0.000005 


0.706789 ± 0.000027 


-968.8 ±0.4 


0.71 


89.42 


A7TH60PH130 


0.950477 ± 0.000004 


0.712722 ±0.000025 


1167.1 ±4.6 


0.71 


129.89 


A7TH60PH210 


0.947876 ± 0.000003 


0.702162 ±0.000029 


2886.7 ±3.0 


0.71 


207.73 


A7TH60PH315 


0.950555 ± 0.000004 


0.712752 ± 0.000030 


-1495.7 ±5.5 


0.71 


316.71 


A7TH120PH0 


0.966806 ± 0.000002 


0.531135 ± 0.000005 


1754.0 ± 6.7 


0.71 





A7TH120PH30 


0.966893 ± 0.000002 


0.531905 ± 0.000005 


1370.2 ±5.4 


0.71 


29.99 


A7TH120PH90 


0.968314 ± 0.000002 


0.537982 ± 0.000005 


-269.0 ± 1.3 


0.71 


86.11 


7TH120PH130 


0.968026 ± 0.000002 


0.536229 ± 0.000005 


-1400.6 ±3.5 


0.71 


129.11 


A7TH120PH210 


0.966893 ± 0.000002 


0.531903 ± 0.000005 


-1370.7 ± 5.4 


0.71 


209.98 


A7TH120PH315 


0.967883 ± 0.000002 


0.535581 ± 0.000005 


1495.8 ±4.1 


0.71 


314.62 


A7TH135PH0 


0.970326 ± 0.000004 


0.494390 ± 0.000006 


1108.0 ± 1.1 


0.71 





A7TH135PH30 


0.969909 ± 0.000004 


0.492329 ± 0.000005 


1328.0 ± 1.9 


0.71 


28.77 


A7TH135PH90 


0.970000 ± 0.000004 


0.493176 ± 0.000006 


775.6 ± 2.0 


0.71 


90.61 


A7TH135PH130 


0.970585 ± 0.000003 


0.495978 ± 0.000006 


-207.1 ±0.6 


0.71 


133.14 


A7TH135PH210 


0.969915 ± 0.000004 


0.492359 ± 0.000006 


-1326.6 ± 1.9 


0.71 


208.32 


A7TH135PH315 


0.970615 ± 0.000004 


0.496080 ± 0.000006 


332.6 ±0.4 


0.71 


318.33 


A9TH15PH0 


0.913073 ± 0.000422 


0.764409 ± 0.004231 


2028.2 ±20.6 


0.90 





A9TH15PH30 


0.912688 ± 0.000396 


0.763674 ± 0.003972 


1764.0 ±23.2 


0.90 


30.06 


A9TH15PH90 


0.913227 ±0.000341 


0.766280 ± 0.003287 


22.9 ± 8.4 


0.90 


91.51 


A9TH15PH130 


0.913683 ± 0.000384 


0.766929 ± 0.003722 


-1323.0 ±9.6 


0.90 


130.98 


A9TH15PH210 


0.912685 ± 0.000394 


0.763651 ± 0.003950 


-1763.9 ±23.2 


0.90 


210.02 


A9TH15PH315 


0.913676 ± 0.000394 


0.766788 ± 0.003835 


1455.4 ± 11.0 


0.90 


316.18 


A9TH30PH0 


0.922385 ± 0.000339 


0.758052 ± 0.005934 


-886.8 ± 13.6 


0.895 





A9TH30PH30 


0.923043 ± 0.000316 


0.760116 ±0.005340 


-2358.8 ±6.0 


0.895 


28.07 


A9TH30PH90 


0.921100 ±0.000166 


0.757236 ± 0.002727 


-3346.5 ±34.1 


0.895 


80.64 


A9TH30PH130 


0.920108 ± 0.000209 


0.753952 ± 0.003666 


-2306.6 ±39.5 


0.895 


122.75 


A9TH30PH210 


0.923043 ± 0.000316 


0.760115 ±0.005341 


2360.4 ±4.9 


0.895 


208.02 


A9TH30PH315 


0.920190 ± 0.000224 


0.753916 ± 0.003973 


2040.8 ±38.7 


0.895 


308.95 


A9TH60PH0 


0.943710 ± 0.000311 


0.740252 ± 0.007605 


3792.7 ±5.4 


0.92 





A9TH60PH30 


0.940226 ± 0.000325 


0.724671 ± 0.005399 


4079.5 ± 10.1 


0.92 


36.83 


A9TH60PH90 


0.945191 ±0.000087 


0.733228 ± 0.001762 


-2352.1 ± 1.6 


0.92 


146.22 


A9TH60PH130 


0.944751 ± 0.000168 


0.735431 ± 0.003341 


-3054.6 ± 0.7 


0.92 


160.22 


A9TH60PH210 


0.940226 ± 0.000325 


0.724671 ± 0.005399 


-4079.5 ± 10.1 


0.92 


216.90 


A9TH60PH315 


0.945098 ± 0.000201 


0.739395 ± 0.004485 


2772.5 ± 0.8 


0.92 


334.02 



provide accurate amplitudes for each velocity compo- 
nent, the above parametrization and fitting to a set of 
full numerical simulations has shown its predictive power 
in a number of occasions; for instance by predicting the 
mass ratio dependence that was later confirmed by sets 
of lengthy numerical simulations [1^. The success of the 



original formula allowed the study of higher-order depen- 
dencies on the spin of the holes. In a previous study [40] . 
we found that the "superkick" recoil (where the two BHs 
have equal mass, equal intrinsic spin magnitudes a, and 
spins lying in the orbital plane in opposite directions) 
has the following dependence on the intrinsic spin a and 




fitted the recoil velocity as a function of the angle 9 that 
the spins make with the orbital angular momentum us- 
ing the data from a series of runs reported in Ref. [53] . 
However, those results are inconclusive since they could 
not model the recoil as a function of <f> or model the pre- 
cession of the orbital plane using the data in Ref. [54] . 

In order to analyze the results of the present simula- 
tions, we use the techniques developed in |46j . Briefly, we 
rotate each configuration such that the trajectories near 
merger overlap. We then calculate the spins in this ro- 
tated frame. The angle <p is then defined to be the angle 
between the AsTHxxPHyyy spin of BH1 (the BH origi- 
nally located on the positive x axis) and the spin of BH1 
in the corresponding AsTHxxPHO configuration. Note 
that, for a given family of fixed spin and spin inclination 
angle 9, the angle (p and 4> differ by a constant, which 
can be absorbed in the fitting constants <\>\ and 4> 3 . We 
then fit the recoil to the 1 form 



FIG. 3: Spin precession frequency versus orbital frequencies 
for a — l/\/2 and a = 0.9 configurations. In each case, only 
the <j>i = configuration is shown. The trend appears to be 
(but see 9 = 30 for a = 0.9) that oj sp i n prcc . increases as 9 
increases, with very little variation with the z component of 
the spin or even total spin a. Note the at late times (larger 
^orbit) i^spin prcc. approaches fi or bit, indicating that the spin 
processes at nearly the same rate as the orbit during the final 
orbit and plunge. 



orientation 4> (the angle between the in-plane spin vector 
and the infall direction near merger), 



V iec = V\ cos(<p - (f>i) + V 3 cos(3ip - 303 



(6) 



V = Vicos(0-0i) + V3Cos(30-3</>3), 
Vi = Vx tl a + V!. 3 a 3 , 

V 3 = V 3A a + V 3 , 3 a 3 , (5) 



where V h3 = (-15.46 ± 2.66) km s" 1 , V 3A = (15.65 ± 
3.01) kin's -1 , and V 3 ,3 = (105.90 ± 4.50) km s"\ while 
Vi i = (3681.77 ± 2.66) km s _1 . From that study, it was 
clear that in the "superkick" configuration, the domi- 
nant contribution, even at large a, is linear in a and 
proportional to cos(</)). Note that because of the small 
contributions of V 3 and Vi t3 , we neglect these terms in 
the statistical studies below (where we take a uniform 
distribution in <f> — (f>\). 

In Ref. [53] an alternative approach to fitting recoil 
and remnant mass and spin of a merged BHB was devel- 
oped. It is based on a Taylor expansion in terms of the 
binary parameters and exploits all the symmetries of the 
problem (note that our approach also incorporates these 
symmetries because it is based on PN formulas for the 
instantaneous recoil) . Using one of our previous set of six 
"superkick" simulations in (TBJ, the authors in [53] fitted 
the recoil velocities to terms in cos(0) and cos(30) to ex- 
tract the cubic (in spin) dependence of the recoil from a 
single set of simulations with constant total spin. We also 
modeled that cubic dependence with more simulations in- 
cluding a range of spin magnitude in [ID] . Ref. [53J also 



for each set of configurations with the same spin and 9, 
and then fit the dominant V\ coefficient as a function 
of (p. Results from these fits are given in Table |IV| and 
Figs. |4j6j Note that A9TH60 runs show the largest dis- 
crepancies in the fit, consistent with the larger errors in 
these simulations due to a coarser global resolution (see 
Sec. Ill above). 

Based on the "superkick" formula d5j), we expected 
that the recoil would have the form 



Vi = (Vi,i + V A a cos 9 + V B a 2 cos 2 9 + V c a 3 cos 3 9) 



a sin 6* 



(7) 



where V\ is the component of the recoil proportional to 
cos0, V\ t i arises from the "superkick" formula, and the 
remaining terms are proportional to linear, quadratic, 
and higher orders in S z /m 2 — a cos 9 (the spin com- 
ponent in the direction of the orbital angular momen- 
tum). Here, we do not consider terms higher-order in 
the in-plane component of A a «2 - qa\ denoted by 
A 1 - (A 1 - oc asin# here), where q is the mass ratio, 
because our previous studies showed that these terms 
were small at 9 = 90°. A fit to this ansatz ^ showed 
that the truncated series appears to converge very slowly 
with coefficients Vi,i = (3677.76 ± 15.17) km s"\ V A = 
(2481.21±67.09) km s" 1 , V B = (1792.45±92.98) km s" 1 , 
V c = (1506.52±286.61) km s" 1 that have relatively large 
uncertainties. In addition, we propose the modification 



Vi 



1 + Ea cos 9 
1 + Fa cos 9 



Da sin 9 



(8) 



which can be thought of as a resummation of Eq. |7]) 
with an additional term Eacos9, and fit to D, E, F 
(where we used the prediction of [40] to model the V\ for 
9 = 90°) and find D = (3684.73 ± 5.67) km s _1 , E = 
0.0705 ±0.0127, and F = -0.6238 ±0.0098. Note that E 
is approximately 1/10 of F, indicating that coefficients in 
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TABLE IV: Fits of recoil velocities as a function of ip for each 
family of configurations with fixed a and 6 to the form Eq. |5j. 



CONF 


Vi(km s" 1 ) 


V 3 (km s- 1 ) 


01 


03 


A7TH22.5 


1764 ± 1 


4.6 ±0.1 


58.36 ±0.01 


280 ± 1 


A7TH45 


2766 ± 1 


41.6 ± 0.9 


203.55 ± 0.02 


83 ±2 


A7TH60 


2972 ± 3 


54 ±3 


342.25 ±0.07 


141 ±4 


A7TH120 


1806 ± 1 


27.4 ± 0.1 


191.78 ±0.01 


62 ± 1 


A7TH135 


1352 ± 1 


16.6 ±0.7 


145.21 ± 0.04 


277 ±4 


A9TH15 


2038 ± 2 


27 ± 2 


178.55 ±0.09 


291 ± 7 


A9TH30 


3408 ± 2 


42 ±2 


284.91 ± 0.03 


285 ±3 


A9TH60 


4171 ± 14 


87 ± 10 


158.3 ±0.7 


17 ± 17 




FIG. 4: A fit of K-ecoil versus <p for the A7TH22 . 5PHyyy (Left) 
and A7TH60PHyyy (Right) configurations. 



this series get progressively smaller faster than in Eq. ([7]). 
Interestingly, a fit to just 



Vi = 



1 



1 + Fa cos 9 



Da sin 9 



failed to produce sensible results (a badly conditioned 
matrix was encountered). The two formulas ^ and Q 
give very similar results for a broad range of a (see Ta- 
bles |VJ|VlJ and|VIl]). We then use Eq. @ to predict the 
recoil for higher spin a = 0.9 and test this formula for 
three angles 9 = 90°, 9 = 60°, and 9 = 15°, with very 
good agreement (see Fig. [6]) . In actuality, both Eq. Q 
and Eq. ([7]) provide accurate predictions for our measured 
recoils at a — 0.9. We show the errors in the predictions 
for the a = 0.9 configurations in Table VI 
We also tried fits to 



1 + Eacos0 + Ga 2 cos 2 9 
1 + Fa cos 9 



Dasinfl, (9) 



but found that the coefficients were not well determined. 
In this case, we found D = 3686.34 ± 8.87, E = 0.055 ± 
0.056, F = -0.638 ± 0.051, G = -0.014 ± 0.050. The 
errors in both E and G in this case are larger than the 
values of the coefficients themselves. We therefore do not 
use Eq. ([9| in the analysis below. Similarly large uncer- 
tainties are encountered if the quadratic F correction is 
in the denominator of Eq. ^ rather than the numerator. 



V. BLACK HOLE SPIN EVOLUTION IN GAS 
RICH GALAXY MERGERS 

Full numerical simulations of BHBs typically start 
when the BHs are at distances of the order of 10M from 




FIG. 5: A fit of Kecoii versus (p for the A9TH15PHyyy (Left) 
and A9TH60PHyyy (Right) configurations. 



5000 



4000 - 



3000 



2000 



1000 




FIG. 6: fit of the recoil (Vi) to the form Eq. Q for the 
a = l/\/2 configurations, and predictions (based on this fit- 
ting) for the a — 0.91 recoils. Note how well the a = 0.91 
curve matches the four measured values. For reference, curves 
corresponding to the original empirical formula prediction 
(which only had terms linear in A) for a = 1/V2 and the 
new formula for a — 1 are also included. Note the skew in 
the velocity profile compared to the linear predictions. 



each other. There are good reasons for this. Numerical 
runs are still extremely expensive, they need to run on 
hundreds of nodes for weeks at a time to obtain accurate 
computations of the gravitational radiation. Those few 
runs allow us to infer generic behaviors of the remnants, 
e.g. the modeling of the recoils by the phenomenologi- 
cal Eq. ([3]) . While these initial separations are extremely 
close by astrophysical standards, most of the nonlinear 
general relativistic effects take place at these and closer 
separations. However, if one wants to study statistical 
distributions of recoils by astrophysical seeds one would 
like to input realistic spin and mass ratio distributions 
for the merging BHB. In a first study of such systems we 
assume an isotropic distribution of the spin direction of 
the BHs [55]. This could represent "dry" binary mergers. 
We point out below the relevance of pre-merger accretion 
to partially align the BH spins with the orbital angular 
momentum. We then perform a preliminary study of an 
extended recoil formula (12 1 to see the differences be- 



tween the predictions of the recoil formula based on only 
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TABLE V: Predictions for a — l/v2 simulations based on 
Eq. (|8| (denoted by pade) and Eq. |7| (denoted by FS), as 
well as the measured V%. Note that the 6 = 90° measured 
value comes from Ref. [40]. Velocities are in units of km s _1 . 



CONF 


Vl, pade 


Vi,FS 




TH22.5 


1760.46 


1754.47 


1764 


TH45 


2771.92 


2777.4 


2766 


TH60 


2967.09 


2967.33 


2972 


TH90* 


2605.5 


2600.57 


2603.4 


TH120 


1802.6 


1811.4 


1806 


TH135 


1354.82 


1348.48 


1352 



TABLE VII: Angle 9 that maximizes the recoil and the max- 
imum recoil as a function of a for Eq. (JsJ) and Eq. |7f. Ve- 
locities are in units of km s~ while angles are measured in 
degrees. 



a 


^pade 


Vpade (km s _1 ) 


6>FS 


Vfs (km s- 1 ) 


0.1 


86° 


369 


86° 


369 


0.5 


70° 


1961 


70° 


1955 


0.707 


62° 


2969 


61° 


2968 


0.91 


54° 


4225 


54° 


4232 


1 


50° 


4926 


51° 


4915 



TABLE VI: Predictions for a = 0.9 simulations based on 
Eq. (|8| (denoted by pade) and Eq. |7| (denoted by FS), as 
well as the measured V%. Note that the 6 — 90° measured 
value comes from Ref. |40) . When applying Eqs. Q and |7|) 
we use Q mC rgcr- The relative error quoted here is the relative 
error in the prediction based on Eq. |S|. Velocities are in 
units of km s _1 . 



CONF 


Vl,pade 


Vl,FS 


Vl ,meas 


Rel. Error 


TH15 


1990 


1905.23 


2038 


-2.4% 


TH30 


3367 


3302.23 


3408 


-1.2% 


TH60 


4251 


4258.62 


4171 


1.9% 


TH90* 


3353.1 


3346.76 


3350.41 


0.1% 



linear terms in the spins and the new updated formula 

m 

The spins of massive BHs binding in binaries as a result 
of galaxy mergers can be deeply affected by gas accretion 
during the last stages of their orbital decay (for separa- 
tions < 100 pc). This is due to the gas overdensities that 
the galaxy mergers are expected to convey into the nuclei 
of the galaxy remnants. Such dense gas structures have 
a disk like morphology ( "circumnuclear disks" ) , reminis- 
cent of the initial net angular momentum of the inflowing 
material, as observed in high resolution simulations [e.g. 
[SHIIST] as well as in real merging systems [e.g. l58rfBU] . 

Bogdanovic et al [24] proposed a physical process that 
could align the BH spins with the angular momentum of 
the nuclear disk in which the binary orbit is embedded, 
thus leading to slow recoils for the BH remnant. The evo- 
lution of the spin directions is due to the torques exerted 
by the gas accreting onto the BHs. Since this process 
happens on a timescale shorter than the orbital decay of 
the BHs in the remnant nucleus, the spins tend to align 
before the two BHs bind in a binary. As a consequence, 
the evolution of the spin of each BH in this earlier phase 
can be studied independently, neglecting the presence of 
the second BH. 

The evolution of spin direction and magnitude is gov- 
erned on small scales (milli-pc, much smaller than the cir- 
cumnuclear disk within which BHs, and their accretion 
disks, are embedded). As shown by 6 lj , if the orbital 
angular momentum of an accretion disk around the BH 
is misaligned with respect to the BH spin, the coupled 



action of viscosity and relativistic Lense-Thirring preces- 
sion ('inertial frame dragging') causes important changes 
in the structure of an accretion disk, warping the disk. 
The equilibrium profile of the warped disk can be com- 
puted by solving the equation: 



1 d . « 1 d ( «dn 
RdR {RL ^-RdR{^ R 



R dR \ 2 dR 



dR 

2G S B n x L 



R? 



(10) 



[see [53], where R is the distance from the BH, vr is the 
radial drift velocity, X is the surface density, il is the 
Keplerian angular velocity of the gas, and v\ (y-i) is the 
radial (vertical) viscosity. L is the local angular momen- 
tum surface density of the disk, defined by its modulus 
L and the unit vector I that defines its direction. The 



disk profile that is a solution of Eq. ( 10 ) is composed of 
three regions, whose relative importance depends on the 
values of the specific disk parameters. In the outermost 
region the angular momentum of the gas is unperturbed 
by any relativistic effect, and therefore the direction of 
the disk's angular momentum is independent of the BH 
spin. In the inner region, the fluid is forced to rotate in 
the equatorial plane of the rotating BH, on either pro- 
grade or retrograde orbits. Therefore in this inner re- 
gion the disk is either completely aligned or completely 
antialigned with respect to the BH spin. Finally, be- 
tween the inner and outer regions, characterized by dif- 
ferent directions of their angular momenta, there exists a 
transition region, centered at ~ 100 — 1000 gravitational 
radii, where the disk is warped connecting the inner and 
outer parts of the disk, misaligned one with respect to 
the other. 

The spin of a BH embedded in a warped disk evolves 
under the influence of the disk itself. The BH spin evo- 
lution is described by the equation |62| : 



dSsH 
dt 



MA(R lso )l(R lso ) 



AttG 



disk 



L x Sbh 
R 2 



dR. 
(11) 

The first term in Eq. (Ill accounts for the angular mo- 
mentum deposited onto the BH by the accreted gas at 
the innermost stable orbit (ISO), where A(i?igo) denotes 
the angular momentum per unit mass [Eq. 12.7.18 in 163] 
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evaluated at i?iso and H-Riso) the local disk angular mo- 
mentum direction, which is parallel to «Sbh as discussed 
above. The second term describes the interaction of the 
BH spin with the warped disk. It is responsible for the 
evolution of the BH spin direction, and tends to align the 
direction of the BH spin with the angular momentum of 
the outer regions of the accretion disk. 

The efficiency of the alignment depends on the dy- 
namics of the inflowing material that fuels the small 
scale accretion disks. It is particularly relevant to de- 
termine whether the accretion flow maintains a nearly 
constant direction of the angular momentum over the 
growth episode (i.e. the accretion is "coherent"), or not. 
Only a substantial amount of gas (1 ~ 10% of the BH 
mass) accreting from the same plane (for both the BHs) 
can significantly align the two spins. In order to con- 
strain the degree of coherency of the gas accreting onto 
the BHs, and to predict the spin configurations in BH 
binaries, Dotti et al [25] performed numerical simula- 
tions of BH pairs in large scale nuclear disks with the 
N-Body/SPH code GADGET [64], upgraded to include 
the accretion physics. Here we give a short summary of 
the initial conditions for the different runs. For a more 
detailed discussion, we defer the reader to Refs. [2511^5] , 

The two BHs are placed in the plane of a massive cir- 
cumnuclear gaseous disk, embedded in a larger stellar 
spheroid. The disk is modeled with w 2 x 10 6 gas par- 
ticles, has a total mass Mnisk = 1O 8 M , and follows a 
Mestel surface density profile £(i?) oc where R is 

the radial distance projected into the disk plane. Dotti 
et al. truncated the disk at an outer radius of 100 pc. 
The massive disk is rotationally supported in R and has 
a vertical thickness of 8 pc. Gas is evolved assuming 
a polytropic equation of state with index 7 = 5/3 or 
7 = 7/5. In the former case, the disk is termed "hot" 
as the temperature is proportional to a higher power of 
density than in the latter class of models ( "cold" cases) . 
The cold case has been shown to provide a good approxi- 
mation to a gas of solar metallicity heated by a starburst 
[66l[67]. The hot case instead corresponds to an adiabatic 
monatomic gas, as if radiative cooling were completely 
suppressed during the merger, for example as a result of 
radiative heating after gas accretion onto the BHs [5o] , 
The spheroidal component (bulge) is modeled with 10 5 
collisionless particles, initially distributed as a P hammer 
sphere with a total mass MB U igc(= 6.98 x Moisk)- The 
mass of the bulge within 100 pc is five times the mass 
of the disk, as suggested by [58]. The BHs are equal in 
mass (tobh = 4 x 10 6 M Q ), and their initial separation 
is 50 pc. A BH is placed at rest at the center of the cir- 
cumnuclear disk, while the other is moving on an initially 
eccentric (eo — 0.7) counter- rotating (retrograde BH) or 
corotating (prograde BH) orbit with respect to the cir- 
cumnuclear disk. Given the large masses of the disk and 
the bulge, the dynamics of the moving BH (secondary) is 
unaffected by the presence of the primary until the BHs 
form a gravitationally bound system. Furthermore, the 
gravitational interaction between the orbiting BH and 



the rotating gas forces the BH to corotate on almost cir- 
cular orbits with the disk [68, 69 , before the BHs bind 
in a binary. As a consequence, the initial orbital con- 
figurations of the BHs do not influence the final degree 
of alignment, that, as will be discussed in the following, 
depends only on the thermodynamical state of the disk. 

To follow the evolution of the BH spins, it is necessary 
to track the dynamics of the gas accreting onto the two 
central objects. In the simulations discussed in [25] a gas 
particle can be accreted by a BH if the following two cri- 
teria are fulfilled: 

• the sum of the kinetic and internal energy of the gas 
particle is lower than 6-times the modulus of its gravita- 
tional energy (all the energies are computed with respect 
to each BH); 

• the total mass accreted per unit time onto the BH every 
timestep is lower than the accretion rate corresponding 
to the Eddington luminosity computed assuming a radia- 
tive efficiency of 10%. 

The parameter b is a constant that defines the degree to 
which a particle is bound to the BH in order to be ac- 
creted. Dotti et al [25] set b = 0.3. Note that due to 
the nature of the above criteria, the gas particles can ac- 
crete onto the BHs only if the time- varying Bondi-Hoylc- 
Lyttleton radius is resolved in the simulations. Such a 
small radius can be resolved only by performing very 
high resolution simulations. The gravitational softening 
parameter of the BHs is 0.1 pc. The gravitational soften- 
ing of the gas particles is set to the same value, in order 
to prevent numerical errors. This is also the spatial reso- 
lution of the hydrodynamical force in the highest density 
regions [84], 

The simulations discussed in Dotti et al [35J cannot 
follow the dynamics of the accreting gas on unresolved 
scales. Dotti and collaborators assume that, below the 
spatial resolution of the runs, gas settles on standard ge- 
ometrically thin/optically thick a-disks [70]. The proper- 
ties of those two unresolved disks (one surrounding each 
of the two BHs of the binary) embedded in the larger 
scale circumnuclear disk, are determined by the proper- 
ties of the accreting material. Each gas particle accreted 
by the BH carries with it mass and angular momentum. 
These are data Dotti et al. used to model the unresolved 
accretion discs around the two BHs, becoming the outer 
boundary conditions for Eq. (10) and (111. Specifically, 



Dotti and collaborators define i e dge as the unit vector 
defining the direction of the angular momentum of the 
accretion disk in its outermost region, i.e. where it is 
unaffected by any general relativistic effect. In a warped 
a-disk the two viscosities (radial, v\ and vertical, 2/2) 
can be described in terms of two different dimensionless 
viscosity parameters, ol\ and through the relations 
^1.2 = cti^Hc s , where H is the disk vertical scale height 
and c s is the sound speed of the gas in the accretion disk. 
Additionally, a 2 = h/i^Oix), with oi\ = 0.1 and / 2 = 0.6 
[71]. Further details on the procedure used to evolve the 
BH spins can be found in [52"] . 

The resulting distributions of spin magnitudes and in- 
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FIG. 7: The probability that the dimensionless spin of a BH 
in a merging binary has a given magnitude a for BHs in cold 
disks (squares) and in hot disks (circles). The fits to Beta 
functions are reasonably good. 




FIG. 8: The probability that the spin of a BH in a merging 
binary is at an inclination angle 9 with respect to the orbital 
angular momentum for cold (squares) and hot (circles) cir- 
cumnuclear disks. The fits to Beta functions are reasonably 
good, but miss the small tail in the cold distribution. 



clinations with respect to the angular momentum of the 
newly formed binary are shown in Figs. [7] and [8j respec- 
tively. Red circles refer to BHs embedded in hot disks, 
blue squares to BHs in cold disks. In both the cases, 
the well defined angular momentum of the large scale 
nuclear disk results in coherent accretion flows onto the 
two BHs, and, as a consequence, in high spins strongly 
aligned with the angular momentum of the BH binary. 
Note that, in absence of any alignment, the distributions 
in Fig. [8] should be oc sin(0) in the whole interval [0,71"]. 
As discussed in |25 , a "hotter" disk, with a stiffer equa- 
tion of state, is more pressure supported in the center, 
and, as a consequence, the degree of alignment is lower. 
Because of this additional support, the accretion rates 
onto the BHs in the hot runs are lower, corresponding to 
spin distributions less skewed towards high spin values. 

In order to perform the statistical studies below, we 
fit the spin magnitude and inclination angle distribu- 
tions above to Beta distributions, which have the form 
P{x) oc (1 — x)^ 1 ) x^ 1 ). Fits to the dimensionless 
spin magnitudes for BHs in hot and cold gaseous envi- 
ronments give a = 3.212 ± 0.258, b = 1.563 ± 0.093, and 
a = 5.935±0.642, b = 1.856±0.146, respectively. A com- 
parison of these fits with the measured probabilities of a 
BH having a given spin magnitude is given in Fig. [7j Fits 
to the inclination angle for the hot and cold cases angular 
distributions give a = 2.018 ± 0.181, b = 5.244 ± 0.604, 
and a = 2.544 ± 0.198, b = 19.527 ± 2.075, respectively. 
Note that these distributions P(9) are for 8 in radians. 
The Beta distribution is not defined for 8 > 1, but the 
data are consistent with near zero probabilities for angles 
larger than 1 radian. A comparison of these fits with the 
measured probabilities of a BH having a given spin di- 
rection is given in Fig. [8] 



VI. EXTENDING THE HANGUP-KICK 
FORMULA 

While accretion will tend to align the spins of the two 
BHs in a BHB with the orbital angular momentum, it 
will not align the in-plane components of the spins. Addi- 
tionally, the expected distribution of mass ratios [72TF74] 
indicates that equal-mass mergers are rare. We therefore 
need a way to extend Eq. Q to generic BHBs. 

Using the same post-Newtonian analysis [ST] as in [50] , 
we can extend formula ^ to less symmetric configura- 
tions by replacing a sin 8 by | a^- — qa^\/(l+q) and a cos 6 
by 2[«2 + ( i ,2 Q ; i]/(l + <z) 2 - Importantly, we are assum- 
ing that terms proportional to |d?^ — qa^] n (for n > 1) 
are negligible. If this is not the case, then our expan- 
sion, which can be thought of as a Fourier sine series, 
would still converge, but our extension would contain er- 
rors that may not be small. For example, if a term like 
(a _L ) 2 a z were present, this would contribute to all even 
components of the Fourier sine series and when extending 
the series, we would have to take this into account. This 
would change the behavior of kick even in more symmet- 
ric configurations. This degeneracy in the interpretation 
of the sine series can be broken by examining configura- 
tions with constant a z (while varying a ) and constant 
a 1 - (while varying a z ). These, and other configurations, 
will be the subject of an upcoming paper. Our justifica- 
tion for not including these terms is that the higher-order 
a 1 - terms are small in the "superkick" configuration. Fur- 
thermore, the accuracy with which formula Q predicts 
the results of our a = 0.91 simulations supports the con- 
clusion that these terms remain small. This can be ver- 
ified by confirming that formula ^ is accurate for all 8 
and a (a subject of our ongoing analysis that will be re- 
ported in a forthcoming paper). We emphasize that the 
proposed extension is an ansatz, that while reasonable as 
a starting point for the modeling, needs to be thoroughly 
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tested and refined. 

Our new ansatz for the recoil velocity modifies Eq. ([3| 
by changing the "superkick" vu term. The ansatz has 
the form (after dropping terms that previous studies in- 
dicated were small [50 J: 

Kecoil ((?)«) = v m e\ +t>j_(cos£ei +sin£e 2 ) +«|| h\\, 
r] 2 (l-q) 



A r 



v± = H 



(1 + 9) 



[1 + B m 77] , 



if 



(i + g) 

U|| = 16/7 2 /(l + g) 



v hl + V A S Z 



V B S Z Z + VcSt 



x \a.2 — qo.\\ cos((/>a — <j>i), 



(12) 



where S = 2(q? 2 + q 2 di)/(l + q) 2 , and the coefficients we 
use in the statistical studies below are H = 6.9 x 10 3 [15] , 
A m = 1.2 X 10 4 , B m = -0.93 12J, and the remaining 
coefficients are obtained from Eq. Q above. 



VII. STATISTICAL STUDIES 



Using Eq. (12 1 and the above fitted spin magnitude 



(a), and direction distributions (9), the mass ratio dis- 
tribution suggested in [72H74] . P(q) cx g~ 03 (l - q), and 
assuming that the two BHs can have arbitrary orienta- 
tions for the in-plane component of the spin (i.e. uniform 
probability in the range < <f> < 2tt), we obtain prob- 
abilities for the recoil velocity magnitude and direction. 
To perform our statistical studies, we choose 10 s configu- 
rations (2 x 10 8 configurations in total) randomly chosen 
based on the above probability distributions and exam- 
ine the predicted recoil magnitude and direction. Our 
results are summarized in Table IVTTT] and Figs. [9] and [TO] 
Although we include the pure unequal mass recoil [v m in 



Eq. (12)], we note that the modeling of the angle £ as a 
function of the binary's parameter is incomplete. How- 
ever v m has normally a nonleading effect. We therefore 
chose a constant £ = 145°, as suggested in our previous 
study [48] . 

In Fig. [TT] we show the probabilities that the recoil has 
a given inclination angle (angle with respect to the or- 
bital angular momentum) for hot and cold disks. Because 
of the 9 —> 180° — 9 symmetry, we map all recoil angles 
to the interval 0° < 9 < 90°. Here P(0) is the proba- 
bility integrating over all possible recoil magnitudes, i.e. 
P{9) = J °° P(9,v)dv. The distribution is normalized 

such that J® P(9)d9 = 1. The angular distribution is 
broader for accretion in cold disks, since that tends to 
suppress the "hangup-kick" and "superkick", while the 
distribution is more sharply peaked near 9 = for hot 
disks. 

This strong angular dependence of the recoil has par- 
ticular relevance for the studies of the observational con- 



TABLE VIII: Recoil velocity probabilities (in percent) for 
BHs in hot and cold disks aligned binaries and the proba- 
bilities for the recoil along the line-of-sight having the given 
magnitude range (denoted by Obs.). For the hot case, there 
is a nontrivial probability of observing a recoil larger than 
2000 km s _1 , but for cold disks, such recoils are suppressed. 
Velocities are in units of km s _1 . 



Vel. (km s' 1 ) 



0-100 

100-200 

200-300 

300-400 

400-500 

500-1000 

1000-1500 

1500-2000 

2000-2500 

2500-3000 

3000-3500 

3500-4000 



(Hot) 



Obs. (Hot) (Cold) Obs. (Cold) 



34.2593 % 
21.1364 % 
11.6901 % 
7.8400 % 
5.7590 % 
14.0283 % 
4.0183 % 
1.0309 % 
0.2047 % 
0.0296 % 
0.0032 % 
0.0002 % 



60.1847 % 
16.9736 % 
8.1110 % 
4.8108 % 
3.0913 % 
5.6593 % 
0.9809 % 
0.1638 % 
0.0223 % 
0.0023 % 
0.0002 % 
4. x 10" 6 % 



41.4482 % 
28.3502 % 
12.503 % 
7.0967 % 
4.2490 % 
5.9309 % 
0.4030 % 
0.0185 % 
0.0005 % 
1 x 10~ 5 % 
0. % 
0.% % 



71.2967 % 
16.8471 % 
6.1508 % 
2.8281 % 
1.3973 % 
1.4258 % 
0.0526 % 
0.0015 % 
2 x ] 
0.% 
0.% 
0.% 
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FIG. 9: Probability distribution P(v) of the recoil magnitude 
for BHBs alignment configurations in hot (top curve, red cir- 
cles) and cold disks BHBs (lower curve, blue squares). The 
velocity is in units of km s _1 . 



sequences of merging and kicked BHs surrounded by pre- 
existing gas disks [THl [75ti80j . The merger of a BHB 
resulting in the remnant BH moving across the matter 
that surrounded the original BHB would greatly affect 
the dynamics of the gas and its thermodynamic state. 
This translates into distinctive electromagnetic signa- 
tures that could reveal the presence of recoiling BHs. 
The effect is very pronounced when the BH recoils in 
the orbital plane with a large magnitude. The strong 
preference for large recoils along the axis of the disk over 
those recoiling along the disk itself can strongly suppress 
the magnitude of such signatures. 

In Fig. [12] we show the recoil velocity distribution for 
hot and cold disks integrated over 15° intervals of the in- 
clination angle 9 (in Figs.[9]and 10 above, the integration 



is over 0° < 9 < 180°). While in Fig. [13] we show the 
integrated probability, Tl(v) of a recoil having velocity 
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FIG. 10: Probability distribution P(v) of the recoil magnitude 
along the line of sight for BHBs in hot (top curve, red circles) 
and cold disks (lower curve, blue squares). The velocity is in 
units of km s _1 . 




FIG. 11: The probability distribution of the inclination angle 
9 of the recoil (measured with respect to the axis of the an- 
gular momentum) for hot (narrower distribution, red circles) 
and cold environments (blue squares). Angles are measured 
in degrees. Note that P(18O°-0) = P(0). These distributions 
were created by mapping 9 — > 180° — 6 for 9 > 90°. 



TABLE IX: Recoil velocity direction (9) probabilities for the 
hot and cold cases. Recoils with 9 > 90° have been remapped 
using the symmetry 9 — > 180° — 9. 



Range 


(Hot) 




(Cold) 


0° - 30° 


61.6% 




47.0% 


30° - 60° 


22.6% 




31.4% 


60° - 90° 


15.7% 




21.6% 


P(v) 


O<0<15 

P(v) 


15<0<3O 






FIG. 12: Recoil velocity distributions for hot (red solid lines) 
and cold (blue dashed lines) dis ks for recoils in angular in- 
tervals 15 degrees wide (Eq. ( 12 l predicts equal probabilities 
for recoils at an angle 9 and 180 — 9). Note how rapidly the 
maximum recoil decreases as a function of 9. Recoils as large 
as 1000 km s _1 must have 9 < 15° and in-plane recoils are 
less than 270 km s - (but see comment about unequal-mass 
recoils in the text). 



v or larger [H(v) = P(v)dv\. If we consider recoils 
within 15° of the orbital axis, we see that velocities up 
to 900 km s _1 are likely (i.e. about 1% probability) for 
cold disks and up to 1600 km s^ 1 for hot disks (see also 
Fig. 13). If we look at angles between 15° and 30° we see 
the recoils are limited to less than 800 km s^ 1 (even for 
hot disks, a recoil of 600 km s _1 has < 0.1% probabil- 
ity). At larger 9 angles, the maximum recoil drops below 
270 km s" 1 . 

Since the most striking effects are likely due to BHs 
recoiling with large magnitudes through the plane of 
the disk, it is interesting to examine the probabilities 
of such events occurring. Even with integrated proba- 
bilities of < 10~ 4 , such events may be observed in large 
surveys of galaxies, e.g. the Sloan Digital Sky Survey 
DR7 contains ~ 9 x 10 5 galaxies if the observable ef- 



fects last long enough. We can see in Figs. [12] and [13] 
that while in hot disks we can observe recoils of nearly 
3000 km s _1 , cold disks limit the maximum observable 
recoil to 2000 km s _1 . On the other end, if we re- 
quire the recoiling hole to be within 30° of the orbital 
plane (i.e. 8 > 60°) we cannot observe recoils larger than 
250 km s" 1 , while at intermediate angles velocities seem 
limited to 400 km s _1 for both cold and hot disks. 

In Fig. |14[ we show the angular distribution for re- 
coils in given velocity ranges. Again, because of the 
9 — > 180° — 9 symmetry, we map all recoil angles to the 
interval < 9 < 90°. For convenience, we plot the prob- 
abilities in degrees rather than radians. For these plots, 
we used 10 8 randomly chosen binaries consistent with 
the above distributions for spin magnitude, spin direc- 
tion, and mass ratio. The maximum angle the recoil can 
make with the orbital angular momentum axis is very 
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FIG. 13: The integrated probability Tl(v) of a recoil having 
velocity v or larger [11(h) = P(v)dv\ for hot (top) and 
cold (bottom) accretion disks for recoils in the ranges 0° < 
9 < 15°, 15° < 9 < 30°, • • • , 75° < 9 < 90°. In both cases, 
for recoils larger than 200 km s - , the recoil probabilities are 
smaller for larger values of 9. In the cold case, low-velocity 
(< 200 km s _1 ) recoils with an angle 0° < 9 < 15° are less 
probable than recoils with angle 15° < 9 < 30°. 



restrictive for large velocities, as is shown in Table |X) 



VIII. DISCUSSION 

We studied in detail a family of BHB configurations 
with full numerical relativity that allowed us to single 
out the "hangup-kick". This effect is expected to be 
relevant in generic BBH mergers since it arises from a 
combination of generic properties of the orbital dynam- 
ics of spinning BHBs, namely the "orbital-hangup ef- 
fect" [81] and the "superkick effect" [IB]. We present 
evidence that this effect increases the maximum recoil 
velocity achievable from the merger of two orbiting BHs 
by up to 1200 km s" 1 with respect to previous estimates, 
approaching nearly 5000 km s _1 . Even more importantly 
this maximum recoil is reached for spins at angles near 
50° with respect to the orbital momentum of the binary 
system. We have also shown evidence that accretion in 
the premerger stage of the binary tends to align spins 
with the angular momentum of the system, leading to 



TABLE X: Maximum recoil angle 9 (angle with respect to 
the orbital angular momentum axis) for given recoil velocity 
ranges. Note here that # max < 5 means that 9 must smaller 
than S or larger than 180° — <5. 



Range 


0max (Hot) 


# m ax (Cold) 


0- 100 km s" 1 


90° 


90° 


100 - 200 km s _1 


90° 


90° 


200 — 300 km s _1 


< 80° 


< 70° 


300 - 400 km s _1 


< 45° 


< 40° 


400 - 500 km s _1 


< 33° 


< 30° 


500 - 600 km s _1 


< 25° 


< 21° 


500 - 1000 km s" 1 


< 25° 


< 21° 


1000 - 1500 km s _1 


< 11° 


< 8° 


1500 - 2000 km s" 1 


< 7° 


< 5° 


2000 - 2500 km s _1 


< 5° 


<4° 


2500 - 3000 km s" 1 


< 4° 


< 2° 


3000 - 3500 km s" 1 


< 3° 






40 60 



FIG. 14: Recoil angle probabilities P(9) for hot (top) and 
cold (bottom) disks aligned binaries. The plots on the left 
show angular probabilities for velocities in the ranges — 
100 km s"\ 100 - 200 km s"\ • • ■ , 500 - 600 km s" 1 . The 
plots on the right show P(9) for velocity ranges of 500 — 
1000 km s~\ • • • , 3000 - 3500 km s" 1 . Probabilities for 9 
and 180° — 9 are equal. In the plots, closed circles correspond 
to the smallest range, followed by squares, diamonds, triangles 
(vertex up), triangles (vertex down), and open circles. The 
circles on the axis are an artifact of the visualization tool. 



distributions that favor the "hangup-kick" configurations 
with respect to the purely in-plane "superkick" ones. Due 
to depletion of nearby matter, the merger itself occurs in 
a "dry" regime where accretion no longer affects the BH 
spins. 

In an attempt to estimate the probability of observ- 
ing such large recoils in real astronomical systems, we 
assumed accretion driven distributions for the spin mag- 
nitudes and directions, based on the two extreme sce- 
narios of cold and hot disks, and assumed a mass ratio 
distribution based on independent estimates, to obtain 
non negligible probabilities of observing recoils of se veral 
thousand km s _1 . In particular, the results in Table VIII 



indicate that surveys exploring peculiar differential radial 
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velocities should observe at least one case of a "line-of- 
sight" velocity above 2000 km s _1 out of four thousand 
merged galaxies (assuming "hot" disks). The probability 
that a remnant BH receives a recoil exceeding the escape 
velocity (in any direction) of giant galaxies (2000 km s _1 ) 
is ten times larger. Probabilities of recoils exceeding the 
escape velocity quickly rise to 5% for galaxies with es- 
cape velocities of 1000 km s _1 and nearly 20% for galax- 
ies with escape velocities of 500 km s _1 . These numbers 
indicate that recoil velocities and modeling the accretion 
of the supermassive BHs in centers of galaxies should 
be important ingredients in understanding the growth of 
supermassive BHs and large scale structure formation in 
the universe. 

Our initial study showed the relevance of recoil and 
accretion modeling in order to better understand how 
BHs evolve and grow in the universe. There are several 
aspects that deserve further study. The recoil formula 
needs to be further tested and developed in the interme- 
diate mass ratio regime and for fully precessing BHBs. 
Accretion needs to be modeled at even smaller scales, 
i.e. at sub-milli-parsecs. In between the accretion regime 
governed by Newtonian physics and the fully nonlinear 
regime, a slow adiabatic inspiral occurs. This interme- 
diate regime can be described by semianalytic methods, 
such as the post-Newtonian approximations. It has been 
pointed out that resonances for certain mass ratios can 
lead to further spin alignment [82, 83J and hence change 



the initial spin distributions used in the applications of 
the recoil formula (12 1. Further theoretical and obser- 



vational explorations of the recoil phenomena are well 
worth being pursued since they could represent the first 
prediction and verification of General Relativity in its 
most highly dynamical and nonlinear regime. 
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